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ABSTRACT 

We investigate the shape of the extinction law in two 1° square fields of the Perseus 
Molecular Cloud complex. We combine deep red-optical (r, ?, and z-band) obser- 
vations obtained using Megacam on the MMT with UKIDSS near-infrared (J, H , 
and if -band) data to measure the colours of background stars. We develop a new 
hierarchical Bayesian statistical model, including measurement error, intrinsic colour 
variation, spectral type, and dust reddening, to simultaneously infer parameters for in- 
dividual stars and characteristics of the population. We implement an efficient MCMC 
algorithm utilising generalised Gibbs sampling to compute coherent probabilistic in- 
ferences. We find a strong correlation between the extinction [Ay) and the slope of 
the extinction law (parameterized by Ry)- Because the majority of the extinction to- 
ward our stars comes from the Perseus molecular cloud, we interpret this correlation 
as evidence of grain growth at moderate optical depths. The extinction law changes 
from the "diffuse" value of Ry 3 to the "dense cloud" value of Ry^ 5 as the column 
density rises from Ay — 2 mags to Ay = 10 mags. This relationship is similar for the 
two regions in our study, despite their different physical conditions, suggesting that 
dust grain growth is a fairly universal process. 

Key words: ISM: dust, extinction - ISM: clouds - methods: statistical 



1 INTRODUCTION 

Interstellar dust affects virtually all astronomical observa- 
tions, yet remains poorly understood. A full prescriptive 
theory which quantitatively explains the production, de- 
struction and steady state conditions of interstella r dust 
does not exist although various simple models (e.g. Ilnoud 
l-Jiil il l have been proposed. In particular, the specific min- 
erals, the shape (or structure), mass/size distribution, and 
total amount of interstellar dust cannot generally be pro- 
duced from first principles. Much progress, though, has 
been made toward a descriptive theory. A basic model 
with separate power-law size distributions of bare graphite 
and silicate grains was shown to adequately reproduce 
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most of the observational evid ence bv lDraine fc Lej l|l984 ). 
IWeingartner fc Draind l|200lf ) updated the model to in- 
clude very small grains — polycyclic aromatic hydrocarbons 
(PAHs) — which are necessary both to reproduce observed 
mid-infrared spectral featur es and to e xplain the transiently 
heated grains identifie d by Sellgrenl (|l984l ). Although the 
IWeingartner fc Draind model is relatively simplis- 

tic (the model uses spheroidal grains although grains are 
probably "fiuffy" and fractal) and slightly ad hoc (the sil- 
icate is "astronomical amorphous silicate" with a few lab- 
measured spectral features modified to accommodate the 
observations) , this model has proven to be quite robust and 
is still widely used to interpret the extinction and emission 
feat ures of dust from the X-ray through the mid-infrared 
fsee lDrainj|2003l . for a review of this model's successes and 
shortcomings) . 

Another nice feature of the IWeingartner fc Draind 



© 2012 RAS 



2 J. Foster, K. Mandel et al. 



l|200lll model is that it directly relates an observational pa- 
rameter to an underlyin g physical change. In particular, the 
highly influential work of lCardelli erall l|l989f) (CCM) found 



that different extinction curves towards different lines of 
sight could all be adequately fit by a one-parameter fam- 
ily of curves: 



Ax = Av{ax + bxR^ 



(1) 



where Rv does double-duty as both the slope of the extinc- 
tion law at V band relating the colour excess E{B — V) to 
the extinction Av 



Rv = 



E{B - V) 



(2) 



and as the parameter distinguishing different extinction 
laws. The coefficients a\ and bx are high-order polyno- 
mial functions of wavelength in the optical, but simple 
pow er laws (with index ~ 1.61) in the n ear-infrared (see 
also iMathij Il 990': Rickc & Lebofsky '1985). The model of 
IWeingartner fc Draine (,2001,) explains this result by observ- 
ing that as grains grow in their model the value of Rv grows 
(in the sense that larger average grain sizes have larger Rv) 
and the extinction law at other wavelengths changes consis- 
tent with the parameterization of CCM. Rv= 3.1 is a typi- 
cal value for the "diffuse" interstellar medium (ISM). Dust 
growth may either be due to the accretion of icy mantles or 
from coagulation of fluffy dust grains. Thus, a connection is 
established between observational changes to the extinction 
law and the physical properties of the underlying model. 

Inside dense clouds, dust undergoes signiflcant pro- 
cessing. The two most important processes are coagulation 
(grains sticking together and growing) and mantle-accretion 
of volatile elements. Some accretion of a volatile mantle onto 
refractory dust grains is presumed to occur in cold dense re- 
gions because observations of gas in molecular cores show 
severe depletion of several important molecules at high den- 
sitie s. This is particularly true for carbon-bearing species 
(see lBergin fc T afalla 2007; . and references therein). Indeed, 
a theory of grains where bare silicate cores are covered by 
carbon-rich mantles c an do a fair job of expl aining the ISM 
dust more generally l|Greenberg fc Lil 119991 ) . Additionally, 
the mid-infrared spectra of dense regions often shows fea- 
tures which are interpreted as ice mantles of volatiles. At 
the same time, coagulation must be important, because the 
extinction per hydrogen atom {Av/N{H)) is observed to de - 
crease in dense regions. This is interpreted bv lMathisI (|l990l ) 
as requiring the coagulation of grains, which can be thought 
of as preventing the interiors of grains from effectively par- 
ticipating in extinction. Adding more material to the extin- 
guishing dust via depletion without also coagulating would 
generally tend to increase the total amount of extinction 
per hydrogen atom simply because there is more material 
blocking radiation. 

Observations of background stars behind a column of 
dust provide the best way to directly measure the shape of 
the extinction law but the wavelength range must be chosen 
with care. Dust exhibits rich spectral features in the mid- 
infrared which are normally ascribed to PAHs. The strength 
of these features is expected to be a function of the pre- 
cise chemical history of a population of dust grains and 
the particular radiation environment to which the dust is 
currently exposed. These features can also produce emis- 



sion which contributes to the mid-infrared fluxes of back- 
ground stars observed through a dust cloud, complicating 
the use of these wavebands and presumably explaining the 
disparate attempts t o nail down this portion of the extinc- 
tion curve (compare Indebetouw et al1l2005l : iFlahertv et al.l 
I2OO7I : [chapman et a l."2009). Wavelengths in the blue-optical 
or UV are quite sensitive probes of the extinction law, but 
it is prohibitively expensive to obtain observations of highly 
reddened stars at these wavelengths. The red-optical and 
near-infrared is thus the best place to study the question. 

A single molecular cloud is a complex hierarchical struc- 
ture, with a dense filamentary web and a vast range of phys- 
ical environments. The cores that eventually form stars are 
often wcU-modcUed as simple Bo nncr-Ebert spheres (e.g. 
Alves et al. 2001; S chnee fc Goodm an 2005). Within such a 
centrally condensed spherical object the observed quantity — 
the total column density, which we shall call Av for con- 
ventional reasons — can be loosely related to the underlying 
physical driver of dust grain growth — the volume density — 
in the sense that increasing column density implies increas- 
ing volume density. Our chain of inference is thus Rv oc (a) 
oc n oc Av, where n is the volume density and (a) represents 
the average size of dust grains. Outside of a single centrally 
condensed object, the connection between column and vol- 
ume density will break down. Nonetheless, because struc- 
tures in molecular clouds tend to be centrally condensed 
this chain of inference may hold, at least for the portions of 
the cloud were we are predominantly seeing through a single 
structure. 

Previous studies of changes in extinction law in molec- 
ular clou ds have mostly b e en at large Av {Av> 10 mag- 
nitudes). ICambresv et al.l (|201ll ) report a transition in 
the extinction law at Av = 20 m agnitudes in Spitzer 
IRAC bands (3.6, 4.5, and 5.8 ^m). |r oman-Ziifiiga et al.l 
( |2007h measured the extinction law at large optical depths 
(up to Av of 60 magnitudes) in B59 and found no evi- 
dence for variation in the extinction law at these depths; 
the favoured extinction law thro ugh the entire cloud was 
the Weingartner fc Drain"^ (|200ll ) model with Rv — 5.5. 
ICambresv et al. I l|2005h report a change in the extinction 
properties around Av = 1 magnitude averaged over the 
galactic anti-centre hemis phere. This an alysis is based on 
a comparison against the iSchlegel et al.l lfl99& ) dust emis- 
sion maps maps, so this could al so correspond to the regime 
where the ISchlegel et al.l ll 19981) dust tem perature correc- 
tion is uncertain ( Arce fc GoodmanI Il999l ) . Studies using 
other methods have also inferred changes in dust proper- 
ties as a function of increasing d epth inside a cloud us ing, 
for example, scattered light (e.g. ISteinacker et al.ll2010l ) or 
changes i n the dust emissivity at submillimetre wavelengths 
(e.g. Ste pnik et al.| [2003). 

In this study, we combine deep multi-wavelength data 
in the red-optical (r, i, z) with near-infrared (J, H, K) data 
from the UKIRT Infrared Deep Sky Survey (UKIDSS) for 
large sections of the Perseus molecular cloud in order to 
test the relationship between Rv and Av over the range of 
column densities where the diffuse ISM-like matter of the 
molecular cloud transitions into the dense-core regime. We 
employ a novel hierarchical Bayesian approach to statisti- 
cally model multiple sources of uncertainty and coherently 
estimate the parameters for individual stars and the popu- 
lation. We find strong evidence that there is a global corre- 
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lation in the sense that the extinction law becomes steeper 
{Rv becomes lar ger which corresponds t o larger dust grains 
in the model of IWeingartner fc Draing l|2001f )) as Av in- 
creases from around 2 magnitudes to about 10 magnitudes 
of extinction. 

In § [5] we describe the acquisition, reduction, and cal- 
ibration of our photometric data. In § [3] we describe our 
model to fit the photometric data using an empirical colour 
locus and a parameterized extinction law. We describe a tra- 
ditional least-squares approach (§[371]) as well as a hierarchi- 
cal Bayesian model (§[3]2]) to overcome some of the problems 
in the least-squares approach. We present the results from 
this hierarchical Bayesian model in § |4l and we conclude in 
§ [5l We describe the mathematical details of our statistical 
model in Appendix A. In Appendix B, we present details of 
our MCMC algorithm for computing statistical inferences. 



2 DATA 

2.1 Megacam Data Acquisition 

Our red-optical data were taken during the second halves of 
the night on November 2nd and 3rd, 2007 using Megacam 
at the 6.5-m MMT on Mt. Hopkinfl We observed two fields 
within the Perseus molecular cloud complex; one field was 
centred on B5 (R.A., Dec = 56.96°, +32.84°; Figure [l]) and 
the other field was cent red on the "West-End" of Perseus 
(see IPineda et al] l2008l . for the definition of this region), 
containing L1448, L1451, and L1455 (R.A., Dec = 51.58°, 
-1-30.50°; Figure[2]). Individual exposures were 80 seconds in 
r, 60 seconds in i and 50 seconds in z. The general pattern 
was to observe all of one field each night, so the West-End 
(WE) was observed on 11/2/07, while B5 was observed on 
11/3/07. However, a few additional WE images were ob- 
tained on 11/3/07 at the end of the night. 

Science observations were interlaced with images of our 
calibration field in order to provide airmass solutions. Our 
calibration field was Per-Cal, a portion of the Perseus super- 
cluster of galaxies included in the Sloan Digital Survey Ex- 
tension f or Galactic Underst anding and Exploration (SDSS- 
SEGUE: lAihara et al.ll201ll ') and centred at R.A., Dec = 
(54.9755°, -1-41.1533°). All r-band science images and cali- 
bration images were taken first during the night, and were 
thus taken between an airmass of 1.0 and 1.06, providing 
insufficient lever arm to fit for an airmass correction term. 
Total exposure times were 800 seconds in r, 600 seconds in i 
and 500 seconds in z. For i and z, data were typically taken 
between 1.1 and 2.0 airmass and the observations were in- 
terleaved with 5 repetitions of the following pattern: i,z,z,i. 
Each field was observed as four different quadrants, with 
some overlap between quadrants (see Figures [T] and [2]) and 
each quadrant was observed contiguously in i and z. Thus, 
for an individual quadrant, the r band images were taken 
near zenith and then the i and z observations were taken in 
an interleaved fashion later in the night. 

The weather was generally good but not perfectly pho- 
tometric. On 11/02/2007 the seeing was excellent with typi- 
cal values of 0.6", sometimes improving to 0.4". Some cirrus 

^ Megacam has since moved to the Magellan Clay telescope at 
Las Campanas Observatory 



clouds were observed during the first half of the night (be- 
fore our data were collected) and were occasionally observed 
during the early part of our observing (mostly during cali- 
bration images of Per-Cal). On 11/03/2007 the seeing was 
relatively poor (0.7 - 1"). Again, a few clouds were observed 
earlier in the night, but then cleared. In B5, the r-band im- 
ages for the first quadrant were taken at 45° with respect to 
north, causing us to have some gaps in our final images and 
catalogue (see Figure [ij. 

2.2 Megacam Data Reduction 

Data reduction of the Megacam data followed a largely stan- 
dard set of pr ocedures in IRAF (Image Reduction and Anal- 
ysis Facility; iTodvl [l993l ) . Bias and fiat frames were made 
and used to process the images. Bad pixels were identified 
and cosmic rays removed. Significant interference fringing 
(where unabsorbed light refiects off the bottom of the CCD 
and interferes with incoming radiation) is present in the i 
and z images, and the de-fringe program by B. McLeod was 
used to remove it. Megacam is made up of 36 individual 
CCDs. Because the de-fringing algorithm scales the inten- 
sity of the fringe pattern before subtracting it from an in- 
dividual CCD, in places where individual bright stars or 
strong surface brightness features dominate the background 
for an entire individual CCD the fringes are poorly sub- 
tracted. Bright stars also cause some "ghosts" or halos on 
the final images, where light scatters internally in the tele- 
scope optics, and these increased surface brightness causes 
additional problems for the de-fringing scaling. Defects from 
bright stars and other features can be seen in both Figures [1] 
&[1 

An illumination correction was also applied, although 
this correction was only significant at r. This correction re- 
moves most large variations in zero-point from one CCD 
chip to the next, but in practice this correction was not suf- 
ficient, and we applied a different zero-point for each CCD 
in the calibration. 

We used SWARP (jBertin et al.ll2002l l to co-add individ- 
ual images (weighting was based on measured noise within 
the images) to produce Figures [1] & (2] and source detec- 
tion was performed on these deep co-added images to use 
in cross-referencing the catalogues. Calibrated photometry 
was performed on each frame individually (so that the air- 
mass and zero-point corrections discussed in the following 
section could be applied) and final source magnitudes were 
the noise-weighted average of the magnitudes measured in 
each frame. 



2.3 Megacam Data Calibration 

Because we wish to fit a parameterized main sequence in 
SDSS colours, it is necessary to place our Megacam pho- 
tometry onto the SDSS (r, i, z) system. This process ac- 
counts for variations in filter and atmospheric transmission 
at different sites, as well as any other optical features of 
the telescope which may change the efficiency with which 
a flux at a given wavelength is measured. Literature colour 
corrections of this type are often based on only a small num- 
ber of comparison s tandard st a rs, an d are thus rather per- 
ilous. For example, ICarpenteil (|200ll) use roughly 50 stars 
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Figu re 1. A three-colour (r, i, z) image of B5 witli contours at 3 and 5 mags of Ay from a 2MASS-based extinction map jRidge et alj 
I2OO6I) . Also shown are the four quadrants used in observing. The most energetic young star, B5-IRS1 is visible at the centre with a small 
reflection nebula around it. Processing artefacts include the "donut" shaped ghosts around bright stars and streaks from satellites. In the 
top-left of the image, green lines show areas not covered by r band due to a mistaken alignment. The bright star in this corner caused 
problems for the de-fringing algorithm, resulting in a corrupt image here. 



to derive the widely-used transformat ion equations between 
Two Micron AH Sky Survey f2MASS: [Skrutskie et al1l2006l ) 
and other near-infrared photometric systems. Our calibra- 
tion fields contain thousands of stars in the SDSS catalogue 
and were taken concurrently with our data, reducing the in- 
fluence of variable atmospheric transparency. This provides 
us with a robust way to calculate these transformations. 

For each filter, x, we solve for CCD chip-dependent pho- 
tometric zero-points (ZP^j) and chip-independent airmass 
(cx) and colour terms (dx) by fitting for 

Inst, = SDSS, -HZP,,j -H (SDSS, -SDSSO (3) 
Inst, = SDSS, -fZP,,j 4- cM + (SDSS, -SDSS,) (4) 
Inst, = SDSS, -HZP,,j -l-c^-f d,(SDSS, - SDSS,) (5) 

where A is the airmass at which a particular image was 
taken. Note that because all our r band images were taken 



at low airmass, we had no leverage on the airmass term in 
this filter and thus do not fit for it. However, since our data 
r band images were also taken at the same low airmass, 
the lack of this term will not significantly influence our cal- 
ibration. We solved for the calibration parameters for each 
night separately to account for the different conditions, and 
list the resulting fits in Table [T] 

The colour terms (dx) proved to be quite significant for 
the r and i bands but insignificant for z. We show the re- 
sults of applying this calibration to the control field stars for 
the second night in Figures |3HS1 for r, i, and z respectively. 
Figures for the first night of data are similar. We experi- 
mented with expanding the colour correction to include a 
term proportional to SDSSi — SDSS, but found that this 
additional complexity did not significantly reduce the ob- 
served scatter around the calibration relation. Plots similar 
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Figure 2. A three-colour (r, ^, 2) image of our d ata from the West-End of Perseus with pink contours at 3 and 5 mags of Ay from a 
2MASS-based extinction map l lRidge et al.ll200(j ). Also shown are the four quadrants used in observing. The three main clouds, L1448, 
L1451, and L1455 are labelled. L1448 and L1455 contain young stars. Visible remaining artefacts include the faint "donut" shaped ghosts 
around bright stars, streaks from some satellites or asteroids, and some residual fringing around the "comet" in the lower-centre of the 
image where the scaling algorithm failed. 



Filter (x) Zero-Points (ZP^j) Airmass (c^;) Colour (d^,) 



WE 



r 


-1.01 - 


-1.08 


N/A 


-0.127 


i 


-0.28 - 


-0.32 


0.0104 


-0.153 


z 


0.8 - 


0.9 


0.073 


-0.0089 






B5 






r 


-1.0 - 


-1.1 


N/A 


-0.133 


i 


-0.21 - 


-0.30 


0.033 


-0.156 


z 


0.9 - 


1.0 


0.093 


-0.0067 



Table 1. Calibration to SDSS. Colour term is SDSSr - SDSS, for 
r, i, z. See Eqns. 3-5 for the definition of these coefficients. 



to Figures I3H5] but versus i-z showed that applying the r — i 
colour correction corrected this colour as well. 

Preliminary results of this analysis were presented by 



[Foster et al.l l|2008l ) , combining these Megacam observations 
with 2MASS data to infer a relationship between Ay and 
Rv similar to what we report here. The photometric uncer- 
tainties of the 2MASS data, however, dominated the error 
budget of this prior analysis and motivated us to modify 
our analysis to incorporate higher quality NIR photometry 
from the UKIDSS survey (see § 12. 4p . The subsequent im- 
provement in our error budget also revealed the presence of 
a photometric offset in the fourth quadrant of the B5 Mega- 
Cam photometry. 

Because the four quadrants for each field overlapped 
slightly, we were able to check the reliability of our calibra- 
tion by comparing the calibrated magnitudes for stars ap- 
pearing in multiple quadrants. The spread in magnitude was 
typically consistent with our estimated measurement and 
calibration uncertainty. The only exception was the fourth 
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Figure 3. The effect of colour correction in r band. Left-panel shows (MMT - SDSS) r-band magnitude after solving for zero-point but 
before correcting for the colour term. Right-panel shows the same difference after applying our colour term. 




Sloon r - I Sloon r - [ 

Figure 4. The effect of colour correction in i band. Left-panel shows (MMT - SDSS) i-band magnitude after solving for zero-point but 
before correcting for the colour term. Right-panel shows the same difference after applying our colour term. 



quadrant in B5. For this quadrant, magnitudes were not con- 
sistent with the other three quadrants (the other three were 
all consistent). We found the median ofTset in the overlap- 
ping stars and offset the magnitudes in the fourth quadrant 
by these amounts (Ar = 0.01 mag, Ai = 0.23 mag, A2: 
= 0.13 mag). Since the i and 2 magnitudes (which were 
observed interleaved together) are most significantly dis- 
crepant, it is likely that some of the cirrus observed earlier in 
the night was in front of our field during these observations. 
We used these offsets with overlapping stars to correct the 
magnitudes in this quadrant. 



2.4 UKIDSS Data 

We obtained UKIDSS Data for our two regions from the 
WFCAM Science Archivfl B5 is covered in the Galactic 
Clusters Survey (GCS) while most of th e West- End is coy - 
ered in the Galactic Plane Survey (GPS: lLucas et al]|2008l ). 
We used the fo llowing quali t y cuts in the SQL query as rec- 
ommended by iLucas et al.l l|2008l ) for high- reliability pho- 
tometry: 

^ |http : //surveys . roe . ac .uk/wsa/ index .html | 
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Figure 5. The effect of colour correction in z band. Left-panel shows (MMT - SDSS) 2-band magnitude after solving for zero-point but 
before correcting for the colour term. Right-panel shows the same difference after applying our colour term. Note that the correction was 
insignificant for this band. 



jAperMag3 > -10 and hAperMag3 > -10 and 
k_lAperMag3 > -10 and jAperMagSErr < 0.03 
and hAperMagSErr < 0.03 and k_lAperMag3Err < 0.03 
and jEll < 0.2 and hEll < 0.2 and k_lEll < 0.2 
and jpperrbits < 256 and hpperrbits < 256 and 
k_lpperrbits < 256 and pstar > 0.99 and 
sqrt(hXi*hXi + hEta*hEta) < 0.3 and 
sqrt(k_lXi*k_lXi + k_lEta*k_lEta) < 0.3 
and mergedClass !=0 and 
(PriOrSec=0 or PriOrSec=f ramesetID) 

and cross-matched the Megacam data with the UKIDSS 
data for stars within 0.6". The UKIDSS survey data is ap- 
proximately on the 2MASS colour system, although the fil- 
ters and system responses are slightly different. We discuss 
the influence this may have on our results in Section 13.41 
The UKIDSS data is preferable to the 2MASS data because 
the uncertainties on 2MASS colours for our stellar sample 
are many times larger than the uncertainties on our Mega- 
cam data and so 2MASS errors would dominate the errors 
in our analysis. 



3 FITTING PROCEDURES 

3.1 Least-Squares Fitting Approach 

The traditional least-squares method for determining cor- 
relation between Ay and Rv involves first deriving best-fit 
estimates for each star's Ay and Rv and then testing these 
best-fit values for correlation with another least-squares pro- 
cedure. 

We formalise the details of our fitting procedure as fol- 
lows. Our input data are five observed colours On derived 
from 6 photometric bands and their associated errors cro„. 
Our set of colours {0„) is (r — i,i — z,z — J, J — H, H — K). 



Our model relies on an empirically-derived, fifth-order poly- 
nomial parameterization of the intrinsic colours of the stellar 
locus, Cn'. 

5 

C„ = ^ Dk,nX' (6) 

where the p a rame ters, -Dfc,„ defining each locus come from 
ICovevet al.1 (|2007l ). and are based on an examination of the 
SDSS and 2MASS colours of 600,0 00 point so urces at low 
extinction values. The equation in lCovev et al. I (|2003) uses 
g — i colour in Eqn. [6]as a proxy for stellar type. We use the 
variable x = g — i instead, to emphasize that we do not use 
g-hand data in this study. Eqn. [S] simply provides a means 
to convert a single parameter (a;, the intrinsic stellar type) 
into our observed colours (the set Cn). 

This stellar locus has finite (observed) width, which 
is due to both intrinsic and measurement dispersion. This 
width is different for each colour, n, and varies along the 
sequence (i.e. is a function of intrinsic stellar type, x). We 
average this dispersion along x to produce a single number 
characterising the width of the observed stellar locus, crc„, 
and we assume that the intrinsic width of this distribution 
is one-half the reported width, to crudely account for the 
difficult-to-assess infiuence of measurement error. 

Our model also includes extinctions. En, where 

En = Ax, - (7) 



Ax, ^ Av{ax,+bx,Rv^) (8) 

and Al represents the effective central wavelength for a par- 
ticular band. The coefficients ax, and bx, are taken to be 
constant for each band, although formally, this is an inte- 
gration across the band which is contingent on the stellar 
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Rv =2.0 
Rv =3.1 
Ry =5.3 
KO star 



Ay =10 mag 



"Ay =1 mag 



Figure 6. Red-optical (r 



ICastelli 



i ptical [r — I versus i — z) colours of a 

Kurucj {gOOl) model KOV star produced in SYNPHOT 



by convolving the stellar spectrum with the filter response curves 
for the MMT. The stellar spectrum is reddened with a CCM Ry 
= 3.1 reddening law in steps of Ay = 1 magnitude (plus symbols). 
The change in the effective wavelength of the reddened spectrum 
causes the colours to trace out a curved line in colour-colour space, 
rather than the straight vectors assumed in our model (straight 
lines). This effect makes it look as if a star behind more extinction 
is being reddened by an extinction law with a smaller Ry . 



spectrum (parameterized in our model by x), the filter re- 
sponse, and the atmospheric opacity at the time {t) of each 
observation : 



A, 



Filter(A) x Atmosphere(A, t) x SED(A,a;)dA. 



(9) 

The extinction in a given waveband is therefore not a sin- 
gle number, but a function of the intrinsic stellar SED and 
the amount of extinction that starlight has already passed 
through. It is mor e correct to think o f reddening curves 
rather than vectors (|Stead fc Hoardl2009l . uses the term 'red- 
dening track'). Figure[6]shows this effect by plotting the r — i 
versus i — z colours of a KO V stellar spectrum with increas- 
ing amounts of extinction assuming a CCM law with Ry = 
3.1 (produced with SYNPHOtQ). The extinction vectors plot- 
ted are appropriate for an un-reddened stellar spectrum, but 
as the spectrum undergoes increasing amounts of extinction 
the effective wavelength changes. As Ay increases, the re- 
sulting stellar colours describe a curve, rather than a straight 
vector. 

Figure [6] shows the worst case in our study, since it con- 
siders the r — i and i — z colours which are most sensitive 
to reddening and includes up to 10 magnitudes of visual ex- 
tinction, which is the maximum value we infer for any stars 
in this study. Most stars are behind only a few magnitudes 
of Ay or less, and the longer-wavelength colours are less in- 
fluenced. Including this effect would increase the complexity 
of the model significantly since we would have to iteratively 
compute the reddening. Therefore we do not include it for 
this study (but see Section [3]4] for a discussion of systematic 
biases and other possible errors), where the column of dust 



[www. stscl . edu/resources/sof tware_hardware/stsdas/synphot] 



Filter 


^^ff (A) 




bA 


T 


6374 


0.930 


-0.240 


i 


7797 


0.799 


-0.533 


z 


9052 


0.675 


-0.619 


J 


12554 


0.414 


-0.380 


H 


16358 


0.260 


-0.239 


K 


22002 


0.162 


-0.149 



Table 2. CCM Coefficients 



is relatively small (< 10 magnitudes of extinction in Ay), 
but it should not be ignored when considering the extinction 
law in denser regions. 

We evaluate ax and b\ from the CCM model. Defining 
^ = X~^fim~^, for ^ > 1.1 (i.e. for r and i) 

y EE e- 1.82, (10) 

where the numerical factor of 1.82 is simply for convenience, 

ax = 1.0-)-0.177y-0.504y^-0.02V-(-0.721y^-f0.020y^-0.775yV0.330/ 

(11) 



bx = 1.413y-f2.283y^-fl.072y^-5.384y''-0.623y^-f5.303y®-2.090/ 



while for ^ < 1.1 (i.e. for z, J, H, K) 



ax = 0.574y 



-0.527?/. 



(12) 
(13) 

(14) 
(15) 



These coefficients are summarised in Table [2] 

We then attempt to minimize (using 
MPFlTFUN(!M arkwardtl |2009| ). aka minpack-1) for each 
star the quantity 

^ 2 

0„^C^{x)^E„{Av,Ry) ' 



M 



E 



(16) 



with the constraints that Ay > -0.5 mags, 2.1 < Ry < 5.5, 
and < X <4.5 (recall that x = g — i and parameterizes 
the spectral type of the star). We allow Ay to go slightly 
negative because the intrinsic Cn are derived from fitting 
real stars with some small amount of Ay. 

The problem with this approach is that we often run 
into limits on Ry. Relaxing the Ry limits does not help; 
there is no local minimum in -space along this dimension. 
In particular, in the case of low Ay, we have very little han- 
dle on Ry, and not a great deal on Ay (the intrinsic spec- 
tral type, X, is normally fit robustly in such a case). The 
large number of points which hit our limits on Ry make 
it hard to test the hypothesis that Cov^Ay ,Ry) 7^ for 
confidence. Another statistical drawback of this approach 
is that applying standard statistical estimators to a sample 
of estimates in the presence of error results in biase d esti- 
mates of the population varia nce and correlation (e.g. iKellvl 
l2007l : lLoredo fc Hendrvl201Gl ). That is, when each individual 
star's [Ay,Ry) parameters are fitted separately, the result- 
ing ensemble of estimates will be wider than the intrinsic 
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Figure 7. Least-squares fits for Ay and Rv for tlie stars in our 
data samples for B5 (left panel) and the West-End (right panel). 
The shades of grey show the density of points at a particular lo- 
cation in the diagram. A representative error bar is shown, which 
displays the median uncertainty on Ay and Ry for a single star 
(rather than for one of the bins displayed in this figure). Stars 
for which no local x^-minimum exists between Ry = 2 and Ry 
= 5.5 "pin" at these extreme values. These points are therefore 
unreliable. 



variance, and the apparent correlation of the ensemble will 
be diminished relative to the true intrinsic correlation. 

We plot these results for both B5 and the West-End in 
Figure [T] Despite the large number of points which hit the 
limit, there is a suggestion of a correlation in the remaining 
points for values of Ay above 1 magnitude. This correlation 
is in the direction expected — larger values of Ry at higher 
column density. 



3.2 A Hierarchical Bayesian Approach 

We developed a hierarchical Bayesian model to address some 
of the problems present in the least-squares fitting approach. 
This model allows us to constrain the properties of the dust 
towards, and the intrinsic colours of, individual stars by si- 
multaneously modelling the populations of dust properties 
and stellar colours. These parameters of the full population 
are called hyper-parameters, and describe the shape of the 
distributions of Ay and Ry. An advantage to this approach 
is that we can include the correlation of Ay and Ry explic- 
itly in our model as a hyper-parameter, and then examine 
its marginal posterior probability distribution conditional 
on our full data-set in order to make probabilistic inferences 
about the dust population. Furthermore, by taking a fully 
Bayesian approach, we can coherently compute the joint un- 
certainties in the estimates of both the hyper-parameters 
governing the population and the dust and colour parame- 
ters of individual stars. By modelling the intrinsic popula- 
tion distribution directly, the hierarchical approach coher- 
ently accounts for the relevant estimation errors when in- 
ferring the intrinsic population covariances, thus correcting 
the biases that plague standard estimators in the presence of 
error. The mathematical details of the statistical model are 
given in AppendixlX] All other aspects of the dust model are 



as described in § 13.11 The statistical and numerical meth- 
ods used here are similar i n principle to those introduced by 
iMandel et alj l|2009l . l201ll ) for estimating dust extinction to 
Type la supernovae using optical and near-infrared photom- 
etry. 

The distribution of column densities {Ay) within a 
molecular cloud is o ften observed to be log-normal (e.g. 
iGoodman et al. 20091). though this is not always true. For 
example. iLombardi et all (|2008l ) fin d a log-normal di s tribu - 
tion for Ophiuchus and Lupus, but iLombardi et al.l (|2006l ) 
see a more complex density distribution for the Pipe neb- 
ula, possibly due to multiple clouds along the line of 
sight. Many simulations of turbulent molecu lar clouds (e.g. 
IVazauez-Semadenilll994l : lOstriker et aLllioOll ) predict a log- 
normal volume density distribution. This is equivalent to a 
log- normal column density distribution if the column density 
along each line of sight through the cloud is dominated by 
a single feature. We therefore take a log-normal distribution 
as a prior on Ay. 



log(ylv-) ~ N{fiA,crA). 



(17) 



The distribution of Ry values within a molecular cloud is 
unknown. However, in the range of column densities probed 
in this study we expect Ry to occupy a range of values 
peaked near the value for the diffuse ISM {Ry — 3.1) or 
slightly larger. A Gaussian in the inverse ofRy {ry = Ry^) 
is a particularly mathematically convenient choice: 



Ty 



N{llr 



(18) 



We assume a fiat (uninformative) prior for the intrinsic spec- 
tral type, X, between 0.2 and 4.2, as our experience suggests 
that this is a parameter well-constrained by the data. See 
§ 13.41 for a further discussion of this choice. Our probabil- 
ity model for the stellar colours is that each measured colour 
comes from a normal distribution, with the mean colour vec- 
tor (conditional on the intrinsic spectral type, x) given by a 
5th order polynomial (Eq. |6}, and the standard deviations 
in each colour given by a 5t h order polynomial derived from 
the uncertainties quoted in lCovev et al.l (|2007l ). 

We further assume the intrinsic spectral type, x, is not 
a priori correlated with the other parameters in the pop- 
ulation, although a weak correlation is possible with Ay 
because only bright blue stars are visible through high col- 
umn density material. We aim to test the hypothesis that 
Ay and ry have a non-zero population correlation, which 
we express as 



ry 
log(^i') 



iV 



fir 
/J-A 



(19) 

In fact our model could be extended to describe a gen- 
eral polynomial dependence between ry and Ay. However, 
we restrict ourselves to a linear dependence, in which case 
the joint population distribution of log{Ay) and ry is suf- 
ficiently described by a linear correlation parameter p. The 
exact form of this linear relationship is 



log Ay = ao + ai{ry - 0.32)/0.04 + e 



(20) 



where ao and ai are the intercept at Ry — 3.1 and the 
slope of the linear relationship, respectively. The numerical 
constants are simply for convenience. The residual variance 
Var[e] = a^, and ex. are related to the mean and covariance 
hyper-parameters of Eq. 1191 via Eqs. IA7I 



© 2012 RAS, MNRAS OOO.mfTSi 



10 J. Foster, K. Mandel et al. 



Dust 
Population 



4* 7?'' 



SpecClass^ 



Intrinsic 
Stellar 
Locus 



T 



AppColors^ 



O, 



■Intr Colors. 



Figure 8. The global posterior density of the unknowns given 
the full data set of apparent stellar colours is represented by a 
directed acyclic graph. Unknown parameters are represented by 
open nodes. Observed data (measured colours O) and knowns 
are represented by shaded nodes. The directed arrows between 
the nodes indicate relations of conditional probability. The hi- 
erarchical model coherently incorporates randomness and uncer- 
tainties due to measurement error (purple arrows), spectral class 
{x = g — i) and the intrinsic stellar locus (/xc(a;), Yici^)) (blue ar- 
rows), and dust extinction and reddening (Ay, Ry)^ and its popu- 
lation {a, (T^, ^r, o"^) (red arrows) into inferences about individual 
stars and the population. The probabilistic graphical model de- 
scribes a conceptual mechanism for generating the observed data. 



Again representing the intrinsic stellar colour as x = 
g — i, and denoting the parameters and data for each of 
A*' stars with the label s, the suite of hyper-parameters as 
H = [/ir, cr^, a, c"^], and the vector of observed colours for 
each star as O^, we seek to compute the global posterior 
distribution 



P({At-,r:,a;.};H|{0=})(x 



n 



X P(H) 



(21) 

where we place uniform (non-informative) priors, -P(H), on 
the hyper-parameters. The derivation of this expression is 
described in detail in Appendix]^ The probabilistic struc- 
ture of the statistical model is shown as a directed acyclic 
graph (DAG) in Fig.O In this graphical model, the unknown 
parameters are depicted with white boxes. The knowns and 
observed data that the inferences are conditioned upon are 
described by grey boxes. The relationships of conditional 
probability between all the variables (specified in detail in 

are shown as directed arrows. The graph can be thought 
of as a generative model, or conceptual mechanism for pro- 
ducing the observed data. Given the intrinsic stellar locus, 
for a given spectral class (xs), a set of intrinsic stellar colours 
is generated. From the dust population, random values of the 
dust parameters {Ay,Ry) are drawn. The intrinsic colours 
and dust effects combine to produce the apparent colours, 
which are then sampled with measurement error to generate 
the observed colours Os. Further information on DAGs and 
examples of t h eir applicat i on in astronom y can be found in 
iMandel et all ()2009l . f2oTl} ): iBishod l|2006h . 

Because we are considering several thousand stars at 
once, there is a large number of parameters and hyper- 
parameters to be estimated jointly {3N + 5). Details of an 
efficient algorithm for solving this problem are presented 
in Appendix [B] A Markov Chain Monte Carlo algorithm 



with generalised Gibbs sampling was used to efficiently ex- 
plore the probability distribution of the model parameters. 
This algorithm was designed to avoid getting stuck at spu- 
rious local maxima with extreme Ry (the problem plaguing 
the least-squares fitting approach), and to negotiate multi- 
ple possible solutions in the global posterior distribution. At 
each step we draw a new value of a given parameter by sam- 
pling from the probability distribution of this parameter, 
conditioned on the current value of all the other parameters 
at this step. In each full cycle we sample the intrinsic colour 
(x), Ay and ry for each star, and then sample the hyper- 
parameters which describe the whole population. After a full 
cycle, the current value of all parameters is recorded as the 
position of the Markov chain. After a sufficient number of 
iterations the Gibbs sampler converges to a stationary solu- 
tion and maps out the full posterior probability distribution 
of parameters and hyper-parameters conditioned on the full 
data set of observed colours. 

We used least square fitting (Eq. I16p to estimate ini- 
tial values for each parameter for each stars. To each value 
we added some random noise to generate different starting 
positions for each individual Markov chain. We remove the 
first 10% of each chain as "burn-in" , during which time the 
chain is still seeking the equilibrium distribution of the pa- 
rameters. We sample 2500 times with 4 parallel, independent 
chains, but record only every 10th sample to reduce the auto- 
correlation within each chain. We measur e the convergence 
of th e model using the Gelman-Rubin ([Gelman fc RubirJ 
Il992h statistic. This compares the variance between chains 
and within chains to estimate the mixing of the chains and 
diagnose their convergence to the equilibrium posterior dis- 
tribution. We consider a chain to have converged if the max- 
imum Gelman-Rubin ratio is less than 1.1. 

This general approach in which we use population infor- 
mation to infer properties of individual stars i n a Bayesian 
contex t is similar to the work reported in iBailer-JonesI 
(j201lf ). Unlike that work, our model does not include paral- 
lax information. We make no attempt to solve for distance, 
considering only colours rather than magnitudes. In addi- 
tion, our model does not consider prior information about 
the expected density of stars of different spectral types; we 
consider all spectral types (g-i) to be equally probable. Our 
study incorporates a more detailed model for the dust ex- 
tinction as that is the focus of this work. 



In some respects thi s work is similar to iKellv et al.l 

(|2012l ). iKellv et al.l (|2012l ) use a hierarchical Bayesian ap- 
proach to infer changes in dust grain properties as a function 
of dust temperature and column density in the Bok Glob- 
ule CB244. Kell y ct al . (2012) examined the dust emissivity 
spectral index (/3) in the far infrared and millimetre and 
found, in contrast to many prior studies, no anti-correlation 
between /3 and temperature, but a strong correlation be- 
tween P and column density. In that work, each independent 
pixel in an emission map was fit individually, and the over- 
all population modelled with hyper-parameters. In our work, 
each individual star is fit individually, and the overall popu- 
lation is modelled with hyper-parameters. Futur e work could 
combine the information about dust emission in iKellv et al.l 
(2012) with the information in this study about dust extinc- 
tion to more tightly constrain models of dust growth. 
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3.3 Testing the Gibbs Sampler 

To test our hierarchical model for sensitivity to factors which 
could lead to spurious correlations, we constructed fake data 
sets and assessed our ability to recover input parameters. 
Un-red dened stars were dr awn from the empirical distribu- 
tion of ICovev et al.1 (|2007l ) using the observed distribution 
of intrinsic g — i colours as the underlying probability dis- 
tribution. This basic g — i colour was used to generate the 
other five colours (r ~ i, i ~ z, z — J , J — H , and H — K) by 
drawing from Gaussian distributions with means and stan- 
dard deviations determined by the stellar colour loci and 
standard deviations tabulated in lCovev et al. 1(12007'). Values 
of \og{Av) and rv were drawn from specified distributions 
with known hyper-parameters and used to redden stellar 
colours according to our COM parameters. Gaussian errors 
of a specified magnitude were then added to simulate the 
observed colours. 

Our primary goal is to constrain p, the correlation be- 
tween \og{Av) and rv- To test this, we generated test data 
sets, each with 2000 stars. Values of \og{Av) and rv were 
drawn from joint normal populations in order to produce 
data with a particular correlation. Because of finite sam- 
pling, these simulations did not produce a population with 
exactly the specified correlation, so we measured the sample 
correlation in these parameters after generation and use that 
value for comparison. The marginal posterior distributions 
on most hyper-parameters were approximately normal, with 
mean values close to (well within the 95% interval) the input 
values. As would be expected, the posterior distributions for 
the parameters of individual stars covered the input values — 
we had to recover the input values of individual stars reason- 
ably well in order to get the population parameters correct. 
Many stars had very broad posterior distributions on Rv, 
as we expected. 

As a basic test of our sensitivity to the mathematical 
form of our prior distributions, we conducted another series 
of tests in which log(Av) and rv were drawn from a bivariate 
uniform distribution with a linear correlation over a sensible 
range [Av from 0.1 to 10 mag, Rv from 2.0 to 5.5). Again, 
p was recovered accurately. 

3.4 Potential Biases 

We know that our model contains some important simplifi- 
cations which may bias our result. In addition, our catalogue 
of stellar colours has been produced using some simplify- 
ing assumptions that could, since we are considering small 
colour shifts, bias our result. We discuss what we believe are 
the most important assumptions and sources of bias here. 

For each filter (r, i ,z ,J ,H ,K) we assume a single ef- 
fective wavelength for all stars when calculating the param- 
eters in Table [21 but this assumption is not precisely true. 
However, the spread in effective wavelengths for the stellar 
population we study (which is dominated by late-type stars) 
is small. The spread in effective wavelength is largest in z, 
where /^Xeff = 15A between an FO and KO star. 

A more significant shift in the effective wavelength is 
produced by reddening itself. As seen in Figure l6l when ex- 
tinction reaches many magnitudes of visual extinction, the 
extinction tracks are actually curves rather than straight 
vectors (as we assume). This introduces a bias in our work. 



but crucially it is the opposite direction of the effect we infer 
in §131 That is, if we were to calculate Rv from Figure l6l for 
the Av ~ 10 point, we would infer a value of Rv < 3.1, closer 
to Rv ~ 2. Since all stellar spectra behave in qualitatively 
the same way under extinction (effective wavelengths be- 
come longer in each band, but more quickly in blue bands), 
this effect would tend to bias us toward inferring that Rv 
decreases with Av, the opposite of what we actually infer. 
The effect is small for the low extinction we study here, and 
thus probably has only a small impact biasing our results 
tow ard finding l ess of a correlation between Rv and Av- 

iNaoi et "all ll2007l ') point out the dangers of photometric 
transformations in studies which attempt to study the ex- 
tinction law (they also find a change in the extinction law 
as a function of incre asing op tical depth). A similar con- 
clusion is reached by iGosling et al.l (|2009l ) in their study 
of near-infrar e d ext inction toward the nuclear bulge and 
I Ken von et al.l (1 19981 ) in their study of Taurus — their results 
are only properly valid in their own photometric system. We 
have attempted to transform our Megacam data onto the 
SDSS system and the UKIDSS data is calibrated onto the 
2MASS system. These transformations are necessary, since 
we use an observed stellar colour locus in SDSS and 2MASS 
colours. Despite the large number of stars we are able to 
use for these transformations, this transformation probably 
remains our largest source of systematic uncertainty. 

Our prior on intrinsic spectral type (parameterized by 
X = g — i) is relatively uninformative. That is, we use a uni- 
form prior between a range of intrinsic g ~ i colours. The 
true distribution of intrinsic spectral types is probably not 
uniform, but we do not have a good estimate for this dis- 
tribution a vnor i- In 13.31 we use th e observed g — i colour 
distribution from lCovev et all (120071 ) when testing the Gibbs 
sampler, but this distribution is for the portion of the Galaxy 
sampled by the low-extinction SDSS footprint; it is not nec- 
essarily appropriate for the particular portion of the Galactic 
stellar population sampled in this study. Even with our un- 
informative prior, the median posterior uncertainty on the 
intrinsic g — i colour is only 0.24, significantly smaller than 
the width of allowed colours in our prior (0.2 < g — i < 4.2). 
We choose to leave the prior relatively uninformative, rather 
than use a potentially incorrect prior. 

We also do not assume that there is an a priori corre- 
lation between intrinsic stellar type and Av- This correla- 
tion may be is present at large column densities where only 
the brightest background stars and foreground stars are de- 
tected. In extinction studies, methods such as NICEST im- 
prove on the NIR colour excess method (NICE) to account 
for the bias that this correlatio n produces in ext inction maps 
(lLombardill2009l '). As noted in lLombardil \200^ . the bias in 
extinction maps only becomes significant at column densi- 
ties where Av > 10 magnitudes, which is the maximum 
extinction we probe. For this reason, we do not attempt to 
model this correlation. 

Additional possible systematic effects include calibra- 
tion problems with our Megacam data as we saw in the 
fourth quadrant of B5 and correlations in the underlying 
population not represented in our model. Our study fits 
slope of the extinction law under the parameterization of the 
CCM model. If the extinction in these portions of Perseus 
is not well described by the CCM model then our fits to the 
CCM model may not be particularly enlightening. A number 
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Figure 9. The mean and standard deviations of the marginal 
posterior distributions on Av and Ky for all the stars in our 
data samples for B5 (left panel) and the West-End (right panel) 
from our hierarchical Bayesian model. The shade of grey shows 
the density of points at a position within the diagram. A represen- 
tative error bar is shown, which displays the median width of the 
posterior probability distribution for Ay and Ky for a single star 
(rather than for one of the bins displayed in this figure) . Note that 
the widths of the posterior probability distributions for Ay and 
Ry vary greatly among individual stars. When we include only 
sources with Ay < 2 magnitudes the correlation between Ay and 
Ry is no longer distinguishable from zero. For these points we are 
unable to infer any underlying correlation and the tight relation- 
ship seen in this diagram is due to the strong correlation inferred 
for stars with Ay > 2 magnitudes. 



of obje cts in our survey m ay not have intrinsic colours lying 
on the lCovev eraP (|2007l ) tracks. These include unresolved 
background galaxies or quasars, embedded young stars and 
brown dwarfs. These objects probably constitute a small 
fractional contamination (~ 1 0%) at the depths studied in 
this work (see the estimates in [Foster et al.ll2008l ). 



4 RESULTS 

We ran our two regions separately, but the results are simi- 
lar. In each case, a few objects failed to converge, producing 
large maximum values of the Gelman- Rubin statistic. Other 
than that, the chains behaved well, with Gelman-Rubin less 
than 1.1 for 99.9% of the objects in both regions. We identi- 
fied the objects which failed to converge, removed them from 
the catalogue, and re-ran the model. Doing this verified that 
the inferences on the hyper-parameters were unaffected by 
these objects. 

The inferences on Ry and Ay are shown in Figure [51 
These plots show the same basic behaviour in both samples, 
with Ry rising with Ay. The points at low Ay are poorly 
constrained by the data from the individual stars, and the 
posterior estimates of Ry and Ay in this regime are strongly 
informed by the population hyper-parameters. 

To assess the reliability of inferences at low Ay we per- 
formed a series of tests in which we progressively removed 
stars with large inferred values of Ay. We used our ini- 
tial catalogue of inferred values for Ay to remove the most 



highly extinguished stars and reran the analysis with these 
subsets of low-extinction stars. The posterior probability dis- 
tribution on p becomes broader until, when considering only 
stars with Ay < 2 magnitudes, the inference on p no longer 
is inconsistent with zero. That is, if we consider only stars 
with Ay < 2 magnitudes, we are not able to claim that a 
significant correlation exists between Ry and Ay. The in- 
ference on Ry is very weak for these stars because there is 
not a significant amount of reddening. Ultimately, therefore, 
the Ry and Ay inferences for these stars is strongly deter- 
mined by the (linear) relationship we assume between ry 
and \og(Ay) and the hyper-parameters which describe this 
relationship. Without the stars at high extinction we would 
have very little information about the relationship between 
Ry and Ay. One should therefore not pay too much at- 
tention to the relationship below Ay of 2 magnitudes. This 
result gives us some confidence that subtle errors in cata- 
logue colours (due to errors in calibration or photometric 
transformations as discussed in § 13. 4|) are not producing a 
spurious relationship; the lack of strong correlation at low 
Ay makes it more likely that the strong correlation we infer 
at high Ay is genuine. 

An alternative approach would be to expand our model 
to allow for a non- linear population relation between log Ay 
and ry, with which the strength of the correlation could 
change as a function of log Av- This would require building 
a more advance statistical model and significant changes to 
the sampling algorithm, which we leave for future work. 

We compute the posterior distributions of individual 
hyper-parameters marginalised over all the other parame- 
ters. These posterior distributions are all roughly normal so 
we summarise them with the median values and standard de- 
viations in Table O In particular, the two regions are found 
to have relatively similar distributions in ry but the West- 
End has a larger average Ay as well as a larger range in 
this parameter. Translating the median values in Table O 
into the conventional units, B5 has a median Ry of 2.6, and 
a median Ay of 2.2 magnitudes while the West-End has a 
median Ry of 2.8 and a median Ay of 4.3 magnitudes. The 
low median in B5 is easily seen in Figure IIOI which shows 
the spatial distribution of inferred Ay values for both re- 
gions. Our solutions for Ay exhibits spatial structure which 
is consistent with known cloud features, as seen in Figure [TOl 

Consistent with the appearance of Figure [UJ the poste- 
rior distribution of p is negative and inconsistent with zero. 
Remember that p is the correlation between ry and log{Ay), 
and so a negative p implies a positive correlation between 
Ay and Ry. The posterior distributions of p are shown in 
Figure fTTI and the results from the two regions are consistent 
with each other. 

The consistency in posterior distributions exists despite 
the differences between the two regi ons. In their c onsid - 
eration of sub-regions within Perseus, iPineda et al.l (|2008l ) 
found that B5 was the most quiescent region with a simple 
density structure and a simple relationship between column 
density and CO intensity, while the West-End had a more 
complicated relationship between column density and CO 
intensity, suggesting significant clumping within the region. 
Later studies have f ound that the quiescent core within B5 
(jPineda et al.l I2OI0I ) is d ominated by a very narrow (5000 
AU wide) filament |pineda ct al. 2011 ). 

IChapman fc Mundvl (|2009l ) proposed that outflows 
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Figure 10. Spatial distribution of inferred Ay in B5 (left panel) and the West-End (right panel) using our hierarchical model with 
extinction coded in colour. Overlain are Ay contours from the COMPLETE NICER map jRidge et al 11200^ ) with contour levels at 
1.6, 2.4, 3.2, 4.0, 4.8, 5.6 and 6.4 magnitudes of Ay. The spatial distribution of stars with high extinction from our Bayesian analysis 
conforms to known cloud structures. The lower density of points in West-End is principally the UKIDSS/GPS data available for this 
region is less deep than the UKIDSS/GCS data for B5. Gaps at the edges can be due to non-uniform coverage either in our Megacam 
observations (e.g. B5) or in the UKIDSS data (e.g. West-End). Gaps in the stellar distribution in regions of high extinction are areas 
where the extinction is too high for background stars to be seen in our Megacam observations. 



# Stars ctt' Ma <^A (^0 cti P 

B5 

3144 0.391±0.002 0.062±0.001 0.34±0.01 0.31±0.01 0.95±0.02 -0.338±0.008 -0.86±0.01 

West-End 

1278 0.363±0.003 0.082±0.002 0.63±0.02 0.38±0.01 0.94±0.02 -0.287±0.008 -0.84±0.01 



Table 3. Median Marginal Posteriors of Hyper-parameters, /ir and Ur are the hyper-parameters describing the normal distribution of 
ry, which is Ry ; fiA and a a describe the normal distribution of log(Av)- cto and ai describe the linear relationship between ry and 
log(Av) as defined in Egn. I20l p is the correlation between ry and log{Ay). 



within two isolated molecular cores may be the cause of re- 
gions where the mid-infrared extinction law is inconsistent 
with dust models as the outflows impact the dust popula- 
tion. As shown in Figure[Tl B5 is dominated by a single dense 
structure, at the heart of which lies B5-IR S1, which drives a 
powerful multi-parsec molecula r outflow (|Ballv et al ] 1 19961 : 
IYu et al.ll999l : lArce et al.l201(t) . Most of the high-extinction 
stars in the B5 data-set come from areas near this outflow. 
The West-End hosts a few impressive outflows in L1448, but 
also a large number of isolated, dense, cores with little or no 
star formation (see Figure [2]). L1448 is mostly outside the 
UKIDSS coverage of this region. Thus the majority of the 
West-End stars at high extinction are found in cores without 
strong outflows. 

Geometrical differences between the two regions could 
also be expected to produce different correlations between 
Ay and Ry. Recall our basic assumption that Ry oc (a) oc 
n oc Ay, where n is the volume density and (a) is the aver- 
age size of dust grains. Increasing the scatter in the connec- 
tion between n and Ay would tend to increase the scatter 
between Ry and Ay and thus decrease any measured cor- 
relation. As described above, B5 is a relatively simple large 
single struct ure. Velocity infor mation from the COMPLETE 
CO survey (|Ridge et al.ll200^ ) and an exhaustive search for 



outfl ows and shells throughout Perseus by l|Arce et al.ll201Cll . 
I2OIII ) suggest that the region is simple in velocity space 
too. Here, the assumption that column density is simply re- 
lated to volume density is probably reasonably robust. In the 
West-End the story is quite different. A large shell appears 
to be disturbing the entire region, and there are multiple 
velocity components toward some positions — and therefore 
likely multiple independent clouds along some lines of sight. 
These additional velocity components have relatively low 
antenna temperatures and so do not contribute much to the 
measured column density along the line of sight, maintaining 
the connection between column density and volume density. 
The consistency of our results for these two regions suggests 
that geometrical effects are not a signiflcant problem. 

We compute a simple estimate of the volume density 
at which we are seeing significant changes in the extinction 
law. By a column density of Ay = 10 magnitudes we have 
typical Ry values of 4-5. lArce et "all (2011) estimates the 
width of B5 to be no more than 0.5 pc. Ten magnitudes 
of visu al extinction corre sponds to roughly 10^^ cm~^ H2 
atoms iBohlin et al.|[l973 ). If we simply assume a constant 
volume density we can calculate a typical volume density at 
10 magnitudes of Ay inside B5 is 6.6 x10^^ cm~^. Becaus e 
of the density substructure within B5 (jPineda et al.llioill ). 
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Figure 11. The marginal posterior distributions for the hyper- 
parameter p, which describes the correlation between ry and 
\og{Av), for both B5 (solid) and the West-End (dashed). The 
two distributions are consistent with each other and inconsistent 
with p = 0. The sense of the correlation is that Ry increases with 
increasing Ay. 



this estimate represents a lower limit on the volume density 
by which point we are seeing clear evidence of grain growth. 



5 CONCLUSION 

We have used r, i, z deep data from Megacam on the MMT 
combined with UKIDSS J, H, K to probe the extinction law 
in the column density range 0.1 mag < Ay < 10 mag using 
a hierarchical Bayesian model. The Bayesian model allows 
us to place a well-motivated prior on Ay (log- normality), 
and work in a consistent framework with the correlation be- 
tween Ay and Ry treated as just another parameter for 
which we compute the posterior distribution. This model is 
implemented using a generalised Gibbs sampler to generate 
Monte Carlo Markov Chains which explore the full joint pa- 
rameter space. The implementation has undergone a series 
of tests using synthetic data similar to our study data, and 
performs well. 

We find evidence that the extinction law changes over a 
relatively narrow range of column densities, rising from an 
Rv ~ 3 at 2 magnitudes of Ay to an Ry ~ 5 by 10 magni- 
tudes of Ay. Below 2 magnitudes of Ay we have insufficient 
sensitivity to infer a relationship between Ry and Ay. The 
two regions in our study, B5 and the West-End lie at oppo- 
site ends of Perseus and have different physical conditions, 
yet both show a strong correlations between Ry and Ay, 
suggesting that the steepening of the extinction curve (most 
likely via grain growth) is a fairly universal process. 

We deliberately do not provide an explicit relation pre- 
dicting the value of Ry for a given value of Ay. These rela- 



tions are slightly different for the two different regions and 
are effectively determined only over a narrow range of Ay. 
In addition, we expect the underlying predictor of Ry to be 
the volume density, n, not Ay, and the relation between n 
and any observed Ay is likely to be highly variable in differ- 
ent regions. This is particularly the case for observations of 
more distant objects, when substantial column density can 
be the result of a long path through diffuse matter. 

Several recent studies have used Sloan Digital Sky 
Survey data to s tudy dust reddening and extinction (e.g. 



Jones et al.ll201ll : ISchlaflv fc Finkbeineill201ll : ISchlaflv et al] 
SDSS data is typically not of sufficient quality in 
high-extinction regions (we examined the portion of Tau- 
rus covered in SDSS-SEGUE) to be useful in constraining 
the extinction law. Future surveys with deeper coverage in 
the red-optical (Pan-STARRs and LSST) and broader sky 
coverage will provide the necessary depth to apply this tech- 
nique to a large number of molecular clouds without the need 
for dedicated deep observations. 
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APPENDIX A: HIERARCHICAL BAYESIAN 
MODEL FOR STELLAR COLOURS AND DUST 

Hierarchical Bayesian analysis is a framework for proba- 
bilistically modelling multiple sources of randomness and 
uncertainty underlying observed data and unifies infer- 
ence for both populations and individuals of those popu- 
lations. Statistical inference with hierarchical models pro- 
vides a principled method of probabilistic de-convolution 
of physically distinct and random effects that are com- 
bined in the observed data. Probabilistic inference allows 
for not only the estimation of each separate effect, but 
also the exploration of the joint uncertainties and trade- 
offs between the multiple effects. It enables the estima- 
tion of the statistical characteristics of an underlying in- 
trinsic population distribution while accounting for the dis- 
tortions in the observed distribution caused by estimation 
error. This correction is called shrinkage. This is accom- 
plished by hierarchical models via partial pooling, which 
combines the individual estimates with population infor- 
mation. Similar issues regarding inferring the intrinsic dis- 
tributions of inferred quantities in the presence of random 
error have been discussed and specific Bayesian techniques 
have been applied bviKellvi t200 7 ): KcUv fc Bechtold (2007[) , 
iHogg. Mvers. fc Bovvl (|201(]| )^ iLoredo fc Hendrvl (|201ol i. 
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iMandel et all (|2009l . I2OIII ). among others, in other astro- 
physical contexts. An excellen t statistical reference for hier- 
archical Bayesian modelling is lGelman et al. I (|2003l ). 

Our statistical model includes a population distribution 
that models the intrinsic stellar locus of colours, a popula- 
tion distribution for the the dust extinction and reddening 
to each star, and incorporates the measurement error for the 
observed colours for each star. Using Bayes' Theorem, prob- 
abilistic estimates for the unknown parameters of individual 
stars, as well as the hyper-parameters of the populations, are 
coherently derived. In the following sections, we build the 
components of our statistical model and then derived the 
global posterior density of the unknown parameters condi- 
tional on the data. 

Al The Observed Colour Likelihood Function 

Let C be a set of linearly independent intrinsic colours of a 
star. For example, C = {r—i,i — z,z — J,J — H,H — Ks). This 
represents the colours of the stars in the absence of dust. 
Let O be the set of observed colours in the same bands 
with estimates of the measurement uncertainty described 
by a covariance matrix So. The dust absorption in a par- 
ticular band F for a given Av and rv = Ry-^ is modelled 
using the CCM law: Ap — Av{ap + hpry), where the coef- 
ficients ap and bp are determined from stellar spectra. Let 
E{Av ,rv) = Av(Aa + Abrv) be a vector of colour excesses 
in the selected colours for a given Ay, ry, and Aa is a fixed 
vector with elements, e.g. (sr — ai), and Ab is defined simi- 
larly. Then, conditional on the intrinsic colours and the dust 
parameters, the observed colours are O — C+E{Av, rv)+e. 
Under the assumption of Gaussian measurement noise, the 
likelihood function is: 

P{0\ C, Av, rv) = N{0\ C + E{Av, rv), So). (Al) 

The multivariate Gaussian probability density in vector y 
with mean fj,y and covariance matrix Sy is denoted by 

N{y\fJ.y,j:y). 

A2 The Intrinsic Colour Population Model 

We can construct a population model for the intrinsic 
colours of stars using the empirical results of Covey et al. 
2007 (C07) based upon analyses of the un-reddened stellar 
locus using SDSS and 2MASS data. They use the intrinsic 
g — i colour as a proxy for spectral class. We shall refer to 
this variable as x = g — i. Conditional on x, we can con- 
struct a normal distribution of colours using the estimated 
means and "pseudo-standard deviations" computed by C07 
for each bin in x. Here x ranges from 0.1 to 4.3. Hence, we 
can write for any individual star 



PiC\x) = N{C\ ficix),-Ec{x)) 



(A2) 



where fj,c{x) and Sc(a;) are known functions of x based 
upon the locus values of Table 1 of C07. Since they do not 
provide estimates of the colour correlations for a given x, 
we set the diagonal elements of the covariance matrix to the 
squared pseudo-standard deviations in each colour, and the 
off-diagonal terms to zero. To complete the intrinsic popu- 
lation model, we should specify a distribution P{x). We find 
it is sufficient to take a conservative approach and assume a 
flat prior on a;: P(3;) oc 1 over the range of x. 



It will be convenient to analytically integrate out the 
intrinsic colours C. The marginal likelihood of the observed 
colours for a single star, given its parameters x, Av and ry 
is then: 

P{0\x,Av,rv) ^ [ dC P{0\C,Av,rv)P{C\x) 



= N[0\ ncix) + E{Av,rv), So + Sc(x)] 

(A3) 

A3 Dust Population Model with Correlations 

If the parameters of dust to a set of stars are drawn from a 
common population, it is advantageous to model that popu- 
lation. We may expect the dust to stars to have some central 
(log Ay ,rv) values with correlated deviations. A reasonable 
choice is a bivariate Gaussian distribution in (log Ay, ry). 



ry 
log(^v) 



iV 



fir 



P <JA(Jr 



2 

(7.4 



This implies lognormal marginal distribution in Ay, 

P{Av\^.A,o-\) = N{logAv\pA,o-\) X Ay'^, 



(A4) 



(A5) 



and a marginal distribution ry ~ N{^r, cr'r). The population 
mean and variance of log Ay are p,A,<y\, respectively. The 
population mean and variance of ry = Ry^ are p,r,<Jr, re- 
spectively. The linear correlation between log Ay and ry is 
p. We additionally limit Ry to lie in the physically plausible 
range between 2 and 5.5 (0.18 <rv < 0.5). 

The above model can also be understood as a Gaussian 
distribution on ry coupled with a mean linear regression of 
log Ay on ry, i.e. ry ~ N{pr,o-f.) and 



logylv-l ry ^ ao + ai{ry - 0.32)/0.04 -I- e 



(A6) 



where the regression coefBcients and residual variance are 
defined by the hyper-parameters. 



Qo = MA + /5— (0.32 — /ir), 



OMp- 



(A7) 



The error term e ~ A'^(0, a ) has variance Var[e] = a = 
a\{l — p^). Using these equations we can work with either 
{a*a, Pi Mr, ""r} Or {qq , Qi , f"^ , Mr , o"r } ^^'^ translate be- 
tween the parameterizations. For implementing the Gibbs 
sampler, we choose the latter parameterization, since it is 
more easily extended to non-linear relations by adding poly- 
nomial terms to Eq. IA6I 



A4 The Hierarchical Posterior Probability 
Density 

To complete the full probability model we must specify the 
hyper-prior density on the hyper-parameters OL,a^,p,r, and 
crj:. We choose standard diffuse non- informative distribu- 
tions: uniform distributions on pr, ol, log cr^ and log cr^. 

We can now construct the hierarchical posterior den- 
sity. Suppose we have Ns stars with observed colours {Os}, 
s = 1 ... TVs . The unknown parameters for each star are 
the intrinsic colours Cs, the spectral type Xs, and the dust 
parameters Ay,ry. We can, however, analytically integrate 
out Cs and use the marginal likelihood, Eq. IA3I Further- 
more, there are several hyper-parameters describing the dust 
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probability model: a, a^, fir and a^. For star s, the poste- 
rior distribution over all remaining parameters conditioning 
on the observed colour data and on the population hyper- 
parameters is 



P(Av,rY,Xs\ CX,a^ , flr,(7r',Os 



<xP{Os\xs,Al,,r'y) X P{xs] 



(A8) 



Since the hyper-parameters are unknown and must be esti- 
mated jointly with the parameters from the data, we must 
compute the joint posterior density over all unknown param- 
eters and hyper-parameters conditioned on the entire data 
set. This is just a product of A'' terms for each star times 
the hyper-prior. 

P{{AY,rv,Xs};a,a'^,flr,(^r \ {Oa}) <x 



Y\ P{Av, rY,Xs\a, cr^./ir, cr^; Os 



X P{a, a , fj,r,(Tr 



(A9) 

All fully Bayesian inferences are based upon the computa- 
tion of this posterior probability distribution. 



APPENDIX B: MCMC SIMULATION OF THE 
POSTERIOR PROBABILITY DENSITY WITH 
GIBBS SAMPLING 

For a data set with A'' ~ 2000 stars, there are 3N + 5 
parameters and hyper-parameters to be estimated jointly 
from the hierarchical posterior distribution, after analyti- 
cally marginalising out the intrinsic colours {Cs}. To do 
this efficiently, we have developed a Markov Chain Monte 
Carlo algorithm based on Gibbs sampling. 

Markov Chain Monte Carlo (MCMC) is a class of al- 
gorithms that generate random draws from an arbitrary 
probability distribution by simulating a stochastic process 
or random walk through parameter space. Gibbs sampling 
is a particular MCMC strategy that uses the information 
in the set of conditional distributions to make the random 
moves in the stochastic process. One parameter at a time is 
updated from its conditional probability density, holding the 
other parameters fixed at their current values. Cycling this 
process through all the parameters repeatedly generates a 
Markov chain that explores parameter space and converges 
to the joint posterior probability density. Once the random 
samples are generated, inferences can be computed using 
these samples. Further theory and technique s of MCMC can 
be fou n d in various textbooks, e. g. iLiul l|2002l ): lGelman" et al.l 
l|2003h : [Robert fc Casellal (|200d ). 

An advantage of using a Gibbs sampling approach by 
sequential sampling of the conditional densities compared to 
a standard Metropolis strategy is that the Gibbs sampling 
approach does not require tuning the jump size of the pro- 
posal distribution. This is important when probing a high- 
dimensional parameter space as we do here; it is practically 
infeasible to optimise the jump size for thousands of param- 
eters. 

A disadvantage of traditional Gibbs sampling is that it 
allows only for orthogonal moves in the parameter space, 
as one parameter is sampled conditional on the fixed cur- 
rent values of the other parameters. This problem can be 



acute if there are strong posterior degeneracies or correla- 
tions between parameters. For example, in this work, we ex- 
pect there to be trade-offs between the intrinsic colours and 
dust reddening in the fit for a single star's observed colours. 
This trade-off will depend upon the shape of the intrinsic 
stellar locus as well as the current estimate of the reddening 
law parameter Ry for the dust to the star. If there exist 
multiple local posterior maxima in parameter space sepa- 
rated along a directions oblique to the parameter axes, the 
orthogonal nature of the Gibbs sampling moves could cause 
the Markov Chain to get stuck at a sub-optimal solution. 

To alleviate these problems, we have included gener- 
alised conditional sampling steps to our MCMC algorithm. 
These steps enable the chain to move along expected degen- 
eracies between parameters that may be oblique with respect 
to the natu ral co-ordinate system defined by the chosen 
parameters (|Liu fc Sabattil I2OO0I : IlIuI |2002|) . This strategy 
allows for non-orthogonal moves through parameter space 
that change several parameters at once. For example, we 
found it advantageous when fitting for each star to simul- 
taneously reduce dust extinction and increase the intrinsic 
colour while adjusting the reddening law. This corresponds 
to an oblique translation through parameter space described 
by Av — Av + y, X — > x + Cx"/, and ry ry + Cry- The coef- 
ficients {cr, Cx), determining the direction of the translation, 
are randomly chosen according to pre-determined probabil- 
ity distributions, and the magnitude of the translation 7 is 
sampled from the posterior density, such that the stationary 
distribution of the chain is left invariant. These generalised 
conditional sampling steps are described below in Step 4. 
Note that the Gibbs sampler without these steps (i.e. using 
only Steps 1-3, 5-7) generates a valid Markov Chain. The 
additional Step 4 is added to speed convergence of the chain 
and help it escape local maxima. 

We start with initial guesses of the unknown parameters 
and hyper-parameters. The current position of the Markov 
Chain is <S = {{Ay ,rY, Xs}; a, , Hr, o-^). Our Gibbs sam- 
pling strategy proceeds as follows. At each step we sample a 
new value of a parameter conditional on all the others, and 
replace the new value in the state vector <S. The first four 
steps are repeated for all stars, conditional on the current 
values of the population hyper-parameters. 



(i) Draw a new Xs from the conditional posterior density 
Pix,\;{Os}) (X P{Os\xs,A'y,r'y) X P{xs). (We use (•) to 
denote all other parameters and hyper-parameters not writ- 
ten explicitly). This is done by evaluating Eg. IA3P on a fine 
grid in Xs and using the inverse-cdf sampling method. 

(ii) Draw a new Ay from the conditional posterior 



P{At.\;{Os}) ^ P{At.\xs,rt.,a,a';Os) 

cc N{Alr\A,VA) X P(Af.|rC.,a,a'), 



where A — VAf3''"'S{x) ^[Os — fj,c{xs)] is the least-squares 
estimate of Ay with variance Va = {P'^'S{x)~^ , and 
(3 ^ Aa + Abry and S(a;) = Sq J^c{xs). The last term 
is given by Eq. IA6I Since this is non-Gaussian, we obtain a 
draw by evaluating this on a fine grid in Ay and using the 
inverse cdf sampling method. 
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(iii) Sample a new from the conditional posterior 

P{rv\-,{Os}) = P{rv\xi,,Av,OL,a-'^,iir,o-l,Oa) 

ocN{r'v\fv,crl)xN{rl,\nr,(jl) (B2) 
X P{A'v\r'v,a,<j^) 

where fv = a^^^Ab^EQ^ (Os - ^J,c{xs) - Af^Aa) 
is the least-squares estimate of rv and = 
(Ab^EQMb(A^)^)"^ is its variance. A new ry is 
generated by evaluating this on a fine grid in ry between 
O.fS and 0.50 and using the inverse-cdf sampling method. 

(iv) Generalised Gibbs sampling. We translate along di- 
rections involving changes in the three parameters: 

Ay Ay + 7, Xs ^ Xs + c^7, ry^rv+ Cr-^. (B3) 

The direction is randomly chosen each time from the follow- 
ing distributions: 

~ iiV(-0.60,0.05^)+iAf(-0.55,0.02^)+iiV(-0.50,0.05^) 



a.\Y ,D,a'^ ^ N{&,V^a^) (Bll) 

(vii) At this point, we have updated every parameter and 
hyper-parameter, and we record the state of the chain S. 
We return to step 1 and iterate the Gibbs sampler with the 
current parameter values many times until convergence. 

We typically run multiple independent chains in parallel 
starting from randomised positi ons and monitor converg ence 
using the Gelman-Rubin ratio (|Gelman fc Rubinlll99a i. 



Cr ~ iA'^(0.03, 0.02^) -f iAf(0.05, 0.02^) 4- iAf(0.08, 0.02^), 

(B5) 

and the sign of c,. is flipped with 50% probability. These dis- 
tributions were chosen via experimentation and were found 
to significantly help the mixing of the chains. The shift 7 is 
sampled from 

P(7) oc P(Av+7: ''v+Cr7, a;s-(-Ca;7| a, ,[ir,o1\Os) (B6) 

by evaluating this density on a fine grid in 7 and using the 
inverse-cdf method. Then with 7,Ca;,Cr chosen, the transla- 
tion Egs. IBSl is performed. 

(v) Steps 1-4 are repeated for all stars. Once all individ- 
ual parameters have been updated, we update the hyper- 
parameters. First, we sample from the joint conditional 
posterior density P(/ir,o-^|-, {Os}) = P{pur\fy ,&!;.^P((j^\s%^, 
where fy is the sample mean of the current ry of all stars, 
and is the sample variance. This is done by drawing a^. 
from a scaled inverse chi-squared distribution and, condi- 
tional on that, drawing \ir from a normal: 

gI\sI ^lviY-^{N -\,8l) (B7) 



[i.r\fv,(yl ^ N{fy,allN). (B8) 

(vi) Next we sample from the joint conditional posterior 
P(Pl,g'^\-,{OsY) = P(pL,a'^\{\o%A^ yy\). This is equiva- 
lent to the Bayesian ordinary linear regression problem. If 
we let Y be vector with elements logy4^, then we have 
Y\oL,a'^ ,{ry^ ~ N{Da,Ia^). The design matrix D has 
rows Ds = [1, (r^ - 0.32)/0.04]. Compute K = (D^£>)"\ 
d = V^D^Y, and 

Sv^jj^^iY-DAfiY^Da). (B9) 

Gibbs sampling proceeds by drawing a new from a scaled 
inverse chi squared distribution, and then, conditional on 
that, sampling new a from a multivariate normal: 

a^lSy ~Inv-x^(iV-2|Sy) (BIO) 
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